library(ConsensusClusterPlus)
library(pheatmap)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

expr 
dim(expr)

 eSet<-  expr[rownames(expr)%in%c(rownames(C1),rownames(C2),rownames(C3) ),]

write.csv(t(eSet),"expression.csv")
 
clinical<-annCol



realdata <- read.csv("expression.csv", row.names = 1)
realdata[1:3,1:6]

library(survival)


Coxoutput=data.frame()
for(i in colnames(realdata[,3:ncol(realdata)])){
  cox <- coxph(Surv(futime, fustat) ~ realdata[,i], data = realdata)
  coxSummary = summary(cox)
  Coxoutput=rbind(Coxoutput,cbind(gene=i,HR=coxSummary$coefficients[,"exp(coef)"],
                                  z=coxSummary$coefficients[,"z"],
                                  pvalue=coxSummary$coefficients[,"Pr(>|z|)"],
                                  lower=coxSummary$conf.int[,3],
                                  upper=coxSummary$conf.int[,4]))
}
for(i in c(2:6)){
  Coxoutput[,i] <- as.numeric(as.vector(Coxoutput[,i]))
}
Coxoutput <- arrange(Coxoutput,pvalue)  %>% #按照p值排序
  filter(pvalue < 0.05) 

#保存到文件
write.csv(Coxoutput,'cox_output.csv', row.names = F)




# 加载用于聚类的目标基因
degs <- Coxoutput$gene


library(ggrisk)
library(survival)
library(rms)

colnames(LIRI)

LIRI<-realdata[,colnames(realdata)%in%c("futime","fustat",Coxoutput$gene)]

fit <- cph(Surv(futime,fustat)~WDR6+MLIP+UFC1+PCBD1+MAMSTR+
             WNK4+PCBD1+PERM1+NICN1+RTKN+FKBP11+JTB+SNAPIN+IFITM3+
             C19orf60+JUP+OBSL1+AZGP1+UBE2I+PTPLA, LIRI)
ggrisk(fit,cutoff.value='cutoff',
       color.A=c(low='#00BFFF',high='#D2691E'),#A图中点的颜色
       color.B=c(code.0='#00BFFF',code.1='#D2691E'), #B图中点的颜色
       color.C=c(low="#00BFFF",median='#FFFAF0',high='#D2691E'))

surv.dat<-LIRI

f = survival::coxph(Surv(futime, fustat) ~ ., data = surv.dat)
#  2.risk point nomgram.points and lp
riskscore = f$linear.predictors

riskscore <- data.frame(riskscore = as.numeric(riskscore),
                        group = ifelse(riskscore >1.03,"1",ifelse(riskscore >-0.1,"2","3")), # 根据中位数分组
                        row.names = rownames(surv.dat),
                        OS.time = realdata$futime,
                        OS = realdata$fustat,
                        stringsAsFactors = F)
Sur <- Surv(riskscore[,3] ,riskscore[,4])
sfit <- survfit(Sur ~ riskscore[,2],data=as.data.frame( riskscore[,2]) )
ggsurvplot(sfit,
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#D14039", "#EA921D","#008ECB"),
           pval.method = T,
           risk.table =T,
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('A','B','C'))+
  labs(x = "Days")


# 提取用于聚类的基因的表达矩阵
comgene <- intersect(rownames(expr), degs) #提取有表达数据的目标基因
comsam <- intersect(colnames(expr), rownames(clinical)) #提取既有表达矩阵又有样本信息的样本
indata <- expr[comgene,comsam] 

indata[1:3,1:3]
dim(indata)
geneclust<-as.numeric(riskscore[,2])

names(geneclust)<-rownames(riskscore)

samorder <- sort(geneclust)


outTab <- NULL
for (i in rownames(indata)) {
  tmp <- as.numeric(indata[i,names(samorder)])
  cor.res <- cor.test(tmp, as.numeric(samorder), method = "pearson")
  outTab <- rbind.data.frame(outTab,
                             data.frame(gene = i,
                                        r = cor.res$estimate,
                                        p = cor.res$p.value,
                                        stringsAsFactors = F),
                             stringsAsFactors = F)
}

# 按相关性正负来分类
outTab$direct <- ifelse(outTab$r > 0, "A","B") # 正相关标为A，否则标为B
outTab <- outTab[order(outTab$r, decreasing = T),]
table(outTab$direct) # 事实上最终用于展示热图的基因数目更少，因为作者进一步采用了随机森林降维

# 把基因与cluster的相关系数r、p值及其所属的基因集（A、B）保存到文件
write.table(outTab,"ouput_ICIsignatureGene.txt",sep = "\t", row.names = F, col.names = T, quote = F)




annCol <- data.frame("Gene cluster" = ifelse(samorder == 1,"A", ifelse(samorder == 2, "B", "C")),
                     row.names = names(samorder),
                     check.names = F,
                     stringsAsFactors = F)
annCol <- cbind.data.frame(annCol,clinical[rownames(annCol),])

annRow <- data.frame("signature gene" = outTab$direct,
                     row.names = outTab$gene,
                     check.names = F,
                     stringsAsFactors = F)


#annCol<-annCol[,-c(2,10)]
#colnames(annCol)[8]<-"ICI_cluster"
#write.csv(annCol,"annCol.csv")
annCol[is.na(annCol) | annCol == ""] <- "N/A"
annColors <- list()
annColors[["Gene cluster"]] <- c("A" = "#D14039", "B" = "#EA921D", "C" = "#008ECB")
annColors[["gender"]] <- c("male"="#79B789", "female" = "#B5262A")
annColors[["vital_status"]] <- c("alive" = "#79B789", "dead" = "#B5262A")
annColors[["age"]] <- colorRampPalette(c("#00FFFF", "#2F4F4F"))(70)
annColors[["stage"]] <- c("Stage II"="#1E90FF","Stage III"="#C71585","Stage IV"="#DAA520")
annColors[["histological_type"]] <- c("epithelioid"="#1E90FF","spindle"="#C71585","mixed"="#DAA520")
annColors[["MetStatus"]] <- c("Metastatic"="#FFFF66","Non-metastatic"="#9C8C2E")
annColors[["chromosome.3.status"]] <- c("disomy"="blue","monosomy"="red")
annColors[["SCNA_Cluster"]]<-c("A"="#D848A9","B"="#D9681F","C"="#8BE36C","D"="#4169E1")
annColors[["GNAQ"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["GNA11"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["EIF1AX"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["SF3B1"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["Clust2"]] <- annColors[["Clust2"]] <- c("Subtype1"="#008ECB","Subtype2"="#EA921D","Subtype3"="#D14039")
annColors[["signature gen"]] <- c("A" = "#929005", "B" = "#179B16B")

annColors


plotdata <- t(scale(t(indata[rownames(annRow), rownames(annCol)])))
plotdata[plotdata > 3] <- 3 # 截断极端值
plotdata[plotdata < -3] <- -3 # 截断极端值

pheatmap(plotdata,
         cluster_rows = F,
         cluster_cols = F,
         show_rownames = F,
         show_colnames = F,
         annotation_row = annRow,
         annotation_col = annCol,
         annotation_colors = annColors,
         color = colorRampPalette(rev(brewer.pal(n = 8, name = "RdBu")))(100))
#dev.copy2pdf(file = "Cluster_gene_Correlation.pdf", width = 10, height = 8) # 保存图像





